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Abstract 

It is very useful and increasingly popular to extract accurate road centerlines from very-high-resolution (VHR) 
remote sensing imagery for various applications, such as road map generation and updating etc. There are three 
shortcomings of current methods: (a) Due to the noise and occlusions (owing to vehicles and trees), most road 
extraction methods bring in heterogeneous classification results; (b) Morphological thinning algorithm is widely used 
to extract road centerlines, while it produces small spurs around the centerlines; (c) Many methods are ineffective to 
extract centerlines around the road intersections. To address the above three issues, we propose a novel method to 
extract smooth and complete road centerlines via three techniques: the multiscale joint collaborative representation 
(MJCR) & graph cuts (GC), tensor voting (TV) & non-maximum suppression (NMS) and fitting based connection 
algorithm. Specifically, a MJCR-GC based road area segmentation method is proposed by incorporating mutiscale 
features and spatial information. In this way, a homogenous road segmentation result is achieved. Then, to obtain 
a smooth and correct road centerline network, a TV-NMS based centerline extraction method is introduced. This 
method not only extracts smooth road centerlines, but also connects the discontinuous road centerlines. Finally, to 
overcome the ineffectiveness of current methods in the road intersection, a fitting based road centerline connection 
algorithm is proposed. As a result, we can get a complete road centerline network. Extensive experiments on 
two datasets demonstrate that our method achieves higher quantitative results, as well as more satisfactory visual 
performances by comparing with state-of-the-art methods. As another contribution, a new and challenging road 
centerline extraction dataset for VHR remote sensing images is made and publicly available for further studies. 
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(TV) & non-maximum suppression (NMS), fitting based connection. 

Guangliang Cheng, Feiyun Zhu, Shiming Xiang and Chunhong Pan are with the National Laboratory of Pattern Reeognition, Institute of 
Automation, Chinese Academy of Sciences (e-mail: {guangliang.cheng, fyzhu, smxiang and chpan}@nlpr.ia.ac.cn). 



2 


I. Introduction 

Accurate road extraction from remote sensing images is an essential preprocessing step for various 
applications, such as vehicle navigation [1], geographic information system (GIS) update [2] and intelligent 
transportation system [3] etc. However, it is time-consuming, costly and tedious to manually label the 
road area from the image. Thus, it is desired to find ways to automatically extract road areas from 
images. Although, recent research on road extraction [4], [5], [6], [7] have been proposed to address this 
challenging task, they are far from mature. 

According to the resolution of sensors, remote sensing images can be classified as low-, median- and 
high-resolution images. Roads on the low- and median-resolution images tend to be lines with small 
width of only one or two pixels. Recent line detection methods, such as dynamic programming and 
snake algorithms [8], template matching [9] and Hough transform [10], have focused on the extraction 
of road networks from these images. Due to recent advances in remote sensors, a large amount of high- 
resolution images become available, which exhibit more details about the earth surface. It is an urgent 
task to extract road from high-resolution images for various applications. However, compared with road 
extraction from low- and median-resolution images, there are a number of difficulties to extract the road 
area from high-resolution remote sensing imagery. First, small objects can be observed and the images 
tend to be affected by noise. Thus the spectral signatures of road become more heterogeneous. Second, 
complex backgrounds and contextual objects, such as trees, buildings and vehicles on the roads, usually 
appear in the high-resolution images. Finally, there are some road-like objects, such as buildings, rivers 
and parking lots, that may be misclassified as roads. 

Most road area extraction methods [11], [12] are based on the pixel-based classification method, which 
brings in heterogeneous results due to the noise and occlusions under vehicles and trees in the VHR remote 
sensing images. For the road centerline extraction methods, many researchers [12], [13], [14] applied 
the morphological thinning algorithm, which is both fast and easy to perform. However, this algorithm 
produces many small spurs, which reduce the correctness and smoothness of the road network. Though 
some regression based centerline extraction methods [15], [16] are introduced to solve this shortcoming, 
they are ineffective to extract the centerlines around the road intersections. 

To overcome the above shortcomings in the existing methods, we propose a novel three-stage based 
method to extract smooth and complete road centerlines from very-high-resolution (VHR) remote sensing 
images: 1) Homogeneous road area segmentation; 2) Smooth and correct road centerline extraction; 3) 
Centerline connection around the road intersections. The proposed method integrates three techniques, 
that is, multiscale joint collaborative representation (MJCR) & graph cuts (GC), tensor voting (TV) & 
non-maximum suppression (NMS) and fitting based connection algorithm. Specifically, to obtain a homo- 
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geneous road area segmentation result, a MJCR-GC based road segmentation algorithm is proposed, which 
incorporates multiscale contextual features and spatial information. Then, a TV-NMS based centerline 
extraction algorithm is put forward to gain smooth and correct centerlines. Finally, to well connect the 
centerlines in the road intersections, a fitting based connection algorithm is introduced. 

The main contributions of the proposed approach are highlighted as follows: 

• MJCR and GC are combined to obtain a homogenous road segmentation result. In MJCR, a novel road 
based feature is firstly proposed, which integrates spectral, structural and contextual road characteristics. 
This feature is in line with the human perception of road recognition. 

• A new TV-NMS based centerline extraction method is introduced to extract the road network. It can not 
only extract smooth and correct road centerlines, but also connect the nearby discontinuous centerlines 
due to unconnected regions in the segmentation result. 

• To overcome the ineffectiveness of the existing road centerline methods in the intersection areas, a 
fitting based connection algorithm is proposed to complete those unconnected centerlines around the 
road intersections. 

• A new and challenging road centerline extraction dataset is publicly available for further studies. It 
contains 30 VHR remote sensing images together with the corresponding centerline reference maps. 
The remainder of this paper is arranged as follows. The related road extraction work is systematically 

reviewed in Section II. In Section III, the details of the proposed road area extraction and centerline 
extraction method are introduced. Experimental evaluations as well as detailed comparisons between our 
method and state-of-the-art methods are provided in Section IV. Finally, the conclusions will be outlined 
in Section V. 


II. Previous Work 

For VHR images, according to the extracted road results, the existing road extraction approaches can 
be classified into two classes: 1) wad area extraction methods, 2) road centerline extraction methods. 

Road area extraction methods mainly depend on image classification and segmentation. Zhang et al. 
[11] proposed an integrated approach that combines k-means, fuzzy logical classifier and shape descriptors 
of angular texture signature. It can separate the roads from the parking lots that have been misclassified 
as roads. A new method for extracting roads based on advanced directional morphological operators was 
proposed in [17], in which Path Openings and Path Closings were introduced to extract structural pixel 
information. Yuan et al. [18] presented an automatic road extraction method for remote sensing images 
based on locally excitatory globally inhibitory oscillator networks. A multistage framework to extract roads 
from the high-resolution multispectral satellite image was introduced by Das et al. [6]. In this method, 
probabilistic support vector machines and salient features were used. 
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Recently, a convolutional neural network based algorithm was introduced to learn features from noisy 
labels in [19], in which the training labels were generated by applying an algorithm trained on a general 
image dataset. Mnih et al. [20] proposed a deep neural network method, which is based on restricted 
boltzmann machines (RBMs), to extract urban road network from high-resolution images. In this method, 
unsupervised pretraining and supervised post-processing were introduced to improve the performance of 
the road detector substantially. A higher-order conditional random field (CRF) model was applied for road 
network extraction by Wegner et al. [21], in which the road prior was represented by higher-order cliques 
that connect sets of superpixels along straight line segments, and the road likelihood was amplified for 
thin chains of superpixels. 

Most popular and successful road centerline extraction methods consist of one or two processing steps: 
1) classification and 2) centerline extraction. Zhu et al. [22] proposed a road centerline extraction approach, 
which is based on the binary-greyscale mathematical morphology and a line segment match algorithm. An 
accurate centerline detection and line width estimation method via radon transform was introduced in [23]. 
Gamba et al. [24] extracted the road centerline with the help of adaptive filtering and perceptual grouping. 
A novel road centerline extraction method was proposed in [12] by integrating mutiscale spectral-structural 
features, support vector machines (SVMs) and morphological thinning algorithm. 

In recent years, Poullis and You [25] proposed a novel vision-based system for automatic road centerline 
extraction. This system integrated Gabor filtering, tensor voting and global optimization using graph- 
cuts into a unified framework. A novel system [26] was introduced to extract road centerline form high 
resolution images, in which probabilistic road center detection, road shape extraction, and graph-theory- 
based road network formation are utilized. Chaudhuri et al. [13] presented a semi-automatic road centerline 
extraction algorithm. To achieve this, directional morphological enhancement and directional segmentation 
were used to extract the road area, then thinning method was applied to extract the road network. An 
automatic road centerline extraction method was introduced by Miao et al. [16], in which potential road 
segments were obtained based on shape features and spectral features, followed by multivariate adaptive 
regression splines to extract road centerlines. 

Shi et al. [14] presented a two-step method for urban road extraction. First, spectral-spatial classification 
and shape features were employed to obtain road segmentation results. Then morphological thinning 
algorithm was used to extract centerlines. An integrated urban main-road centerline detection method was 
introduced in [15]. Road extraction result was achieved by fusing spectral-spatial classification and local 
Geary’s C method. Then, to extract smooth road centerlines, local linear kernel smoothing regression 
algorithm was introduced. It relieves the shortcoming of thinning algorithm, while it can’t work well 
in the intersection areas. Hu et al. [27] proposed a three-step road centerline extraction method. First, 
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Stagel: MJCR-GC based Road Area Segmentation 



Multiscale over-segmentation 


Road likelihood 




Road segments 



Shape-based 
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Final segment 



Stage2: TV-NMS based Road Centerline Extraction 


Fig. 1. Flowchart of our method. It contains three stages: road area segmentation, road centerline extraction, and road centerline connection. 
In the stage, multiscale joint collaborative representation (MJCR) and graph cuts (GC) are applied to obtain the road segmentation result, 
then a shape based elimination algorithm is introduced to distinguish the road segments and road-like segments. For the 2^^ stage, tensor 
voting (TV), non-maximum suppression (NMS) are proposed to extract smooth and correct road centerlines. A fitting based connection 
algorithm is utilized to complete the road centerline network in the 3^^ stage. 


adaptive mean shift was used to detect road center points. Then, the salient linear features were enhanced 
via stick tensor voting. Finally, a weighted Hough transform was applied to extract the arc primitives of 
the road centerlines. Sironi et al. [28] provided a new idea by setting the road centerline extraction task 
as a regression problem, which learns a scale space distance transform from the high-resolution image. 

III. Road Segmentation and Centerline Extraction 

In this section, we propose a novel urban road centerline extraction method, which consists of three 
stages: road area segmentation, road centerline extraction and road centerline connection. The flow chart 
is shown in Fig. 1. Speciflcally, in the 1®' stage, MJCR and GC are combined to gain an accurate road 
area segmentation result. Then, to eliminate those road-like segments from the segmentation result, a 
novel shape based elimination algorithm is introduced. To obtain a smooth and correct road centerline 
network, the 2"*^ stage provides the tensor voting (TV) and non-maximum suppression (NMS) based road 
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centerline extraction method. In the 3’’^’ stage, to obtain a complete road centerline network, a fitting based 


connection algorithm is utilized around the road intersections. 

A. MJCR-GC based Road Area Segmentation 

In this subsection, we propose two techniques, i.e. multiscale joint collaborative representation and 
graph cuts (MJCR-GC), to classify the image into road and non-road areas. Specifically, there are three 
stages. First, we present a novel object-based spectral, structural and contextual (SSC) method to extract 
features. Second, a MJCR-GC based framework is proposed to obtain an accurate road area segmentation 
result. Third, to remove the non-road objects (i.e. houses and parking lots) from the road class, a new 
shape-based elimination algorithm will be employed. 

1 ) Object-based SSC feature extraction: In the VHR remote sensing image, the road segments tend 
to be elongated and under the occlusions of cars and trees etc. Based on this observation, we employ 
an object-oriented SSC feature to achieve the following two goals: (1) to reduce the side influence of 
occlusions and the spectral variations; (2) to extract the geometric characteristics of road segments. The 
basic assumption behind the object-based algorithm is that spatially adjacent pixels are prone to be grouped 
into spectrally homogeneous objects. Specifically, the simple linear iterative clustering (SLIC) [29] is 
employed to obtain the multiscale over-segmentation for VHR remote sensing images. SLIC adapts a 
local k-means clustering approach according to their color similarity and proximity in the image plane. 
There are three advantages to apply the SLIC for the over-segmentation task. SLIC is not only more 
efficient, but also achieving more satisfactory over-segmentation results compared with the state-of-the-art 
methods, for example graph-based algorithms [30], [31] and gradient-ascent based algorithms [32], [33] 
etc. Besides, there is only one parameter specifies the number of superpixels, which makes it extremely 
easy to use. 

We obtain multiscale segmentations with different superpixel numbers. Intuitively, superpixel number 
controls the segment size. A large superpixel number will averagely result in a small object size, and vice 
versa. As Fig. 1 shows, we get three different over-segmentation results for each image. 

For each over-segmentation result, we define the spectral attribute of an object as the average spectral 
value within this object. For the Lth object, the spectral attribute (SA) [12] is given by 



( 1 ) 



where p is a pixel in the Lth object (1 < p < Ni) and SP’"(p) is the spectral value of pixel p in band m; 
The superscripts r, g and b represent the RGB channels, respectively. To exploit the structural attributes 
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Fig. 2. The illustration of how to extract the spectral, structural and contextual feature for the superpixel with dark border. 


for each object, we utilize the shape index (SI) [12], [34] and aspect ratio (AR) as follows 

Per, 


SL = 


AR. = 


4\/Areai ’ 
Length- 


( 2 ) 


(3) 


Widthi ’ 

where Per^ and Area, denote the perimeter and area of the Ath object, respectively; Lengthy and Width, 
are the length and width of the minimum bounding box surrounding the z-th object. 

The SI measures the smoothness of the object border and the AR describes the length-width ratio. 
Intuitively road regions tend to be elongated ribbons with large perimeters and small areas. Therefore, the 
SI and AR of road regions tend to be large. 

The aforementioned spectral and structural attributes for each object can be combined as a hybrid 
feature (HF), which is defined as 

HF* = {SAf,SI^ARf}, (4) 


where HF| is the hybrid feature for the z-th object at the s-th scale. 

To enhance the discriminative power of each object, the spatially adjacent objects should be considered. 
In the proposed method, co-occurrence among objects are used as the high-level contextual feature. 
Intuitively, there are at least two road segments (i.e. left and right, or up and down) in the neighborhood 
of one road segment. This is the prior information for humans to recognize the road areas under the 
occlusions of trees and shadows. Our SSC feature employs this recognition scheme. Fig. 2 illustrates how 
to extract the SSC feature for one center segment. There are three steps. First, we find out the neighboring 
segments (i.e. the blue border regions) for the center segment (i.e. the dark border region). Second, we 
rank the neighboring segments according to the their similarity values to the center segment, where we 










use the Euclidean distance as the similarity measure. Third, the HFs of both center segment and its top 
two neighboring segments are stacked as the contextual feature for the center segment. 

2) MJCR-GC based road area segmentation: In this subsection, MJCR and GC are combined to 
obtain an accurate road area segmentation result. Specifically, SSC feature is utilized by MJCR to gain the 
road probability of each segment. Then, to enhance the label consistency among the neighboring pixels, 
GC is introduced to incorporate the spatial information. 

a) MJCR based road likelihood: Recently, sparse representation classification (SRC) [35] has 
been proposed for face recognition. SRC represents a testing sample by a sparse linear combination of 
training samples with £i-norm constraint. In the remote sensing imagery, Chen et al. [36] applied a sparse 
framework for the hyperspectral image classification. A similar approach to SRC is the collaborative 
representation classification (CRC) [37], [38]. CRC also represents a testing sample with the linear 
combination of training samples. However, contrary to the f'l-norm regularization in SRC, CRC employs 
an (' 2 -norm regularization. It provides competitive accuracy while with significantly lower computational 
complexity. 

Consider a dataset with n training samples X = where Xj e and d is the feature 

dimensionality. Let y* G {1,2, ...C} be the class label, where C is the number of classes. For the road 
extraction task, we set (7 = 2. Let ni and n 2 be the number of training samples for road class and 
non-road class, respectively. 

For a testing sample x, the corresponding combination coefficient a, based on all the training samples, 
could be obtained via 

a* = argmin ||x - Xq;|| 2 -f A ||rx,sa:|| 2 , (5) 

OL 


where Tx.x is a biasing Tikhonov matrix between the test sample and all the training samples; A is a 
global regularization parameter that balances the representation loss and the regularization term. Note that 
a* is the optimal representation vector of a with n x 1 elements. Specifically, Tx.x e is designed 
in the following form: 

||x-Xi||2 ••• 0 


r 


X,x — 


( 6 ) 


0 ||x-x„||2 

Tx.x is a diagonal matrix, whose diag value measures the discrepancy between a certain training sample and 
the testing sample. Intuitively, if the testing sample belongs to the road class, the discrepancies between the 
testing sample and those road-class training samples are small, while inversely the discrepancies are large 
for the non-road-class training samples. Given a large regularization parameter A, to achieve a minimum 
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objective in Eq. (5), the road-class testing sample are more likely to be represented by road-class samples 
rather than non-road-class samples, thus a* tend to be sparse. 

The representation coefficient a* can be estimated in a closed-form solution as 

a* = (X^X + AT^ ~rx,s) ~'x'^X. (7) 

After that, the training samples X are partitioned as road-class samples Xi and non-road samples 
X 2 , and the coefficient vector a* are partitioned as and accordingly. The residual between the 
approximation and the testing sample can be defined as 


Ri{x.) = llXiO:;* - x| 

ie{i,2} 


( 8 ) 


In this research, we get the class likelihood rather than the class label of each sample. Therefore, we 
define the road-class likelihood as 


Pr(x) = 


-R2(x) 


(9) 


i?i(x) + i? 2 (x)’ 
and the non-road-class likelihood as Pnr(x) = 1 — Pr(x). 

After obtaining the road likelihood for each object at a certain scale, all the pixels in the same object 
are given an identity likelihood value as the object. Then, we fuse all road likelihood maps from the three 
different scales into an integrated one. Finally, the road likelihood of each pixel Xi can be defined as 


Pr{xi) = max p^^{xi). (10) 

se{l,2,3} 

We also tried other fusion strategies, such as mean, median and min. Their results are inferior to the max 
rule. Thus we use the max-fusion rule in our experiments. 

b) GC based road area segmentation: In the remote sensing image, some road areas are under 
conditions of spectral variability as well as occlusions of trees and cars. To relieve the side effect of these 
conditions and get a coherence road extraction result, graph cuts (GC) algorithm is used for the road 
extraction task. 

Given an image /, the GC algorithm constructs an undirected graph Q = {V, £}, where V denotes the 
pixel set in the image and £ represents the set of undirected graph edges between neighbouring pixels 
[39]. For the road extraction task, we define the label “1” for the road class and “0” for the non-road 
class. GC tries to minimize the following objective 

C{C) = Cr{C)^aC,{C), (11) 

where £ is a labeling set, Cr{C) and ChiC) denote the regional term and boundary term, respectively; 
a is a trade-off parameter balancing the two terms. In the road extraction problem, the regional term 



10 



Fig. 3. Comparison of elimination results between our ellipse-based algorithm and the bounding box based algorithm, (a) Original image, 
(b) Final segmentation result, (c) Elimination result via the bounding box based algorithm, (d) Elimination result by ellipse-based algorithm, 
(e) and (f) are the close ups of the rectangle region in (d) and (c) respectively, (g) The illustration of the ellipse-based algorithm and the 
bounding box based algorithm. 


Cr{C) defines the individual penalty to classify each pixel into the road class. The boundary term 
describes the coherence between spatially neighboring pixels. 

We have obtained the road likelihood map via aforementioned MJCR. The regional term is defined as 


Cr{C) = ^ - \Og{pr{Xi)) , ( 12 ) 

i&I 

where Pr{xi) is the road likelihood probability of pixel Xi. Intuitively, the spatially adjacent pixels tend 
to belong to the same class, thus the boundary term is defined to measure the label discontinuity among 
neighboring pixels, which is defined as 


C,{L)= Y, 



(13) 


where J\f denotes a standard 8-neighborhood system, which contains all the unordered pairs of neigh¬ 
bouring pixels. L^.) is the distance metric between the label and L^., i.e. if Lx^ and L^. have 

different labels, we denote = 1, otherwise, we define it as 0. Xi and Xj are the RGB feature 

vectors of pixel Xi and Xj. 11-112 denotes the ^ 2 -norm. To avoid a zero divisor, we add a small value e 
(typically e = 0.001) to the denominator. 

For binary labeling problem, the objective function in Eq. (11) can achieve the optimal solution via the 
mincut/maxflow algorithm in polynomial time [40]. As Fig. 1 shows, coherence road segmentation result 
can be obtained after the GC algorithm. 

As Fig. 1 shows, coherence road segmentation result can be obtained after the GC algorithm. However, 
it remains some road-like objects. Some strategies should be taken to remove those road-like objects. 






(a) 


(b) 


Fig. 4. Illustration of tensor voting algorithm, (a) voting process. The voter (V) casts its votes to the receiver (V). (b) The sketch map of 
different saliency values corresponding to different eigenvectors for a tensor. 


3) Shape based elimination algorithm: After the process of the GC based road area segmentation, 
there still remains some road-like objects (see Fig. 3(b)). To address this issue, the road shape features 
should be employed to distinguish the potential road segments and road-like segments. 

In general, roads have the following characteristics: 1 ) Roads are connected and do not have small 
areas; In our experiments, we regard all the segments whose pixel number is less than a threshold K as 
non-road class. 2) Roads are elongated ribbons with long length and short width. Linear feature index 
(LFI) is defined to evaluate this characteristic, which can be denoted as 


LFI = 


■^box 


w, 


box 


(14) 


where Lbox and Wbox are the length and width of the bounding box (the red rectangle in Fig. 3(g)), 
respectively. Intuitively, road segments have large LFI values. However, the bounding box based algorithm 
may be failed as shown in Figs. 3(c) and 3(f), some road segments (the sketch map in Fig. 3(g)) may 
have small LFI values. 

To overcome the shortcoming above, we propose a new ellipse-based road elimination algorithm (see 
the green ellipse in Fig. 3(g)). We use an ellipse to match each remaining segment after the area constraint. 
We define the new LFI as 


LFIe = 



(15) 


where Le and We are the major axis and minor axis of ellipse, respectively. In our experiments, we set 
the LFIe as 3. Comparing Figs. 3(c), 3(d), 3(e) and 3(f), the real road segment is removed in Figs. 3(c) 
and 3(f) via bounding box constraint, while our ellipse-based elimination algorithm keeps it remained in 
Figs. 3(d) and 3(e). Thus, it demonstrates that our proposed algorithm is more effective than the bounding 
box method. 
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B. TV-NMS based Road Centerline Extraction 


After the above processes, there are still two issues with the gained road segmentation results. Due to 
the image noise and occlusions caused by trees and cars, the extracted road area always have an unsatisfied 
property, that is, discontinuities. Besides, for some real problems, we only consider whether there exists 
a road, and do not care about the width of the road. To relieve the side effect of discontinuity and to 
be more suitable for application, road centerline extraction has been an active research. The commonly 
used road centerline extraction algorithm is morphological thinning, which is fast and easy to implement. 
However, thinning based road centerline extraction algorithm always produces many spurs and bring in 
many false positives, which reduce the smoothness and correctness of road network. To solve this problem, 
TV-NMS is introduced to extract road centerlines, which can produce smooth centerlines and complete 
the discontinuous road centerlines. 

Tensor voting (TV) was originally proposed by Guy and Medioni [41], and later presented a lot of 
improved versions [42], [43]. It is a technique for perceptual grouping and robust extraction of lines 
and curves from images. TV aims at making the local feature measurements more robust by taking the 
neighboring measurements into consideration. In 2-D tensor voting, local image feature are encoded into 
a tensor filed, in which each tensor S at a certain location can be denoted as a symmetric, positive 
semi-definite matrix representation as 


S = 


^xx ^xy 


— AiOiOi^ -|- X 2 ^ 2^‘2 ) 


(16) 


Sxy Syy 

where Ai and A 2 are the nonnegative eigenvalues (Ai > A 2 > 0 ) of the matrix S; ei and 62 are their 
corresponding eigenvectors. As Fig. 4(b) shows, the graphical illustration of the tensor S is an ellipse. 
The tensor S can be decomposed as 


S = (Ai 


^2)6161^ -|- A2 



+ 6262 



(17) 


where eiOi^ corresponds to a degenerate elongated ellipsoid, termed the stick tensor, (eiCi^ -f 6262 ^) 
corresponds to a circular disk, termed the ball tensor, (Ai — A 2 ) and A 2 are the saliency values of the 
stick tensor and ball tensor respectively. The spectrum theorem [44] states that any 2-D tensor can be 
decomposed into the linear combination of stick tensor and ball tensor. 

In the voting process, as Fig. 4(a) shows, assume point V is the voter and point R is the receiver, OO 
is the osculating circle passing V and R. The second order vote is a stick tensor and has a normal lying 
along the radius of the osculating circle at R [45]. According to the Gestalt principles, the magnitude of 
the vote should be a function of proximity and smooth continuation. Thus the saliency decay function 
can be defined as 


DF{s, K, a) = e 


( 18 ) 
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NMS 


(d) 


connection 

(e) 

Fig. 5. The flowchart of the proposed centerline extraction method, (a) Road segmentation result, (b) Stick saliency map after tensor voting, 
(c) Junction map after tensor voting, (d) Incomplete centerline map after non-maximum suppression (NMS). (e) Complete centerline map 
after the fltting based centerline connection. 

where s = is the arc length of VR; I is the distance between the voter and receiver; 0 is the angle 
between the tangent of the osculating circle at the voter and the line going through the voter and receiver 
(see Fig. 4(a)); k = is the curvature, and c controls the degree of the decay with curvature; a is the 
scale of voting, which determines the effective neighborhood size. In the experiment, the parameter c is 
a function of the scale [45] as 

, -16(ff - 1) X log(O.l) 

TT^ 

In this research, the voting process includes two kinds of voting: sparse voting and dense voting. 
The sparse voting describes a pass of voting where votes are cast to locations that contain tokens only. 
The dense voting denotes a pass of voting from the tokens to all the locations within the neighborhood 
regardless of the presence of the tokens. Then the receivers accumulate the votes cast to them by tensor 
addition. 

After tensor voting, according to the saliency value, the point feature can be classified into three types 
(i.e. curve point, junction point and outlier) [41] as follows, 

curve point if (Ai — A 2 ) > A 2 

junction point if Ai « A 2 > 0 • (20) 

outlier if Ai ~ 0, A 2 ~ 0 

v 

According to the above rules, we get the junction map (see Fig. 5(c)) and stick saliency map (see Fig. 5(b)). 

One great advantage of TV is that it only has one parameter (i.e. scale factor a) to be determined by 
users. Actually, in most cases, the scale factor can be set according to the road width, which will be 
discussed later. 
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(a) (b) 


Fig. 6. Illustration of the proposed fitting based connection method, (a) The sketch map of detecting the middle points, (b) The final map 
after connection (the purple lines). 


To get the consistent road centerline, Non-Maximum Suppression (NMS) is applied on the resulting 
stick saliency map (see Fig. 5(b)). The NMS keeps only the locations that correspond to a local maximum 
along a line perpendicular to the local orientation within a neighborhood of width. As Fig. 5 shows, the 
TV-NMS method can extract smooth road centerline, while it is ineffective around the road intersection. 
To tackle this issue, some other measurements should be taken to complete the road centerline network. 

C. Fitting based Road Centerline Connection 

TV-NMS based centerline extraction method can extract single-pixel-width road centerline well. How¬ 
ever, it can’t extract the centerlines around the road intersections properly (see the initial map in Fig. 1). 
To overcome this shortcoming, a new fitting based centerline connection method is proposed. 

We get the junction map after the TV algorithm. As Fig. 6(a) shows, for each junction point (the red 
circle), we search its local area within a certain distance defined by the radius R. For example, there 
are three line segments in the radius areas in Fig. 6(a). For each segment in the radius, we calculate its 
middle location of all the centerline pixels, termed middle point. Then, we link the middle point and the 
corresponding center point of junction area (see the purple lines in Fig. 6(b)). In this subsection, we employ 
the linear fitting algorithm to connect the discontinuous road centerlines around the intersection area. It 
is behind the assumption that road is straight at the intersection area, which is satisfied in most cases. As 
Fig. 5(e) shows, our fitting based method can link the centerlines well around the road intersection area. 

IV. Experiments and Evaluation 

In this paper, the definition of the “VHR” refers to the image with spatial resolution of 1-2 m per pixel. 
The road width in VHR image is 8-15 pixels. The corresponding road reference map is manually labeled. 

To verify the effectiveness of the proposed method, extensive experiments, on the road centerline 
extraction from VHR remote sensing images, have been conducted on two datasets. The proposed method 
is also compared with other state-of-the-art methods. 
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A. Datasets description 

This section introduces the information of two VHR image datasets for the road extraction experiments. 
It should be noted that few VHR urban road datasets are publicly available. Thus we collected 30 VHR 
images from Google Earth and labeled the road reference map and centerline reference map by a hand 
drawing method. This dataset will be publicly available for research. 

Data #1: This dataset contains 30 VHR urban road images with a spatial resolution of 1.2m per pixel. 
The road width is about 12-15 pixels. There are at least 800 x 800 pixels in each image. Most of the 
images are under the conditions of complex backgrounds and occlusions due to trees and cars. In addition, 
there are road-like objects (i.e. houses and parking lots) in the images. All these factors make it very 
challenging to obtain a satisfying road extraction results. 

Data #2: This dataset is a publicly available road dataset' provided by E. Turetken et al. [46]. It 
contains 14 images with a spatial resolution of 1.5m per pixel. The road width in this dataset is about 
8-10 pixels. Some images are under the conditions of complex backgrounds and occlusions of trees. We 
manually labeled the centerline reference map of each image in this dataset. 

B. Compared algorithms 

To verify the performance, the proposed method is compared with four related methods. The main 
information of all these methods are summarized as follows: 

1) Huang’s method (Huang): Huang et al. [12] introduce a novel road detection method based on 
multiscale structural features and support vector machine. Then, morphological thinning algorithm is used 
to extract the road centerlines. 

2) Miao’s method (Miao): Miao et al. [16] present a road extraction method based on spectral and 
shape features. Then, to overcome the shortcomings of morphological thinning algorithm, multivariate 
adaptive regression splines is introduced to get the smooth and correct road centerlines. 

3) Shi’s method a (Shi“): This method [14] integrates the spectral-spatial classification, shape features 
and local Geary’s C to extract road network. Then a morphological thinning algorithm is applied to extract 
centerlines. 

4) Shi’s method b (Shi'’): To get the road network, Shi et al. [15] fuse the spectral-spatial feature and 
homogeneous property via support vector machine. After that, a local linear kernel smoothing regression 
algorithm is introduced to extract the centerline. It is a state-of-the-art method for road centerline extraction. 

' http://cvlab.epfl.ch/data/delin 
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Fig. 7. Visual comparisons of road area extraction results with four comparing methods. There are two rows and six columns of subfigures. 
Each row shows results on one dataset. From left to right: (a) original image, (b) result of Huang [12], (c) result of Miao [16], (d) result of 
Shi“ [14], (e) result of Shi^ [15], (f) result of Proposed^. From (b) to (f), true positive (TP) is marked in green, false positive (FP) in red, 
false negative (FN) in blue. (Best viewed in color) 


5) Proposed method with three scales (Proposed^): As Fig. 1 shows, the image is oversegmented 
with three different numbers of superpixels, such as 8000, 10000 and 12000. Then the following steps in 
Fig. 1 are performed to get the road centerline. 

6) Proposed method with one scale (Proposed^: To investigate the effectiveness of multi-scale 
information fusion, the image is oversegmented with only one certain number of superpixels, for example 
10000. Other steps are performed as Proposed^ states. 

It should be noted that the codes of the other four comparing methods are not readily available in the 
Internet. We implement them according to the details described in their papers. Our implemented codes 
can achieve the same performance as the original paper presents. In the following experiments, we adjust 
the parameters to gain the best performance of the four methods for fair comparison. Among the four 
comparing methods, regression based centerline extraction algorithm is used in Miao and Shf’, thus they 
can obtain smooth centerline network. Meanwhile, Huang and 5/?/“ utilize the morphological thinning 
algorithm to extract road centerlines, which will produce some small spurs. 

C. Evaluation metrics 

To comprehensively evaluate the proposed method, comparison experiments are conducted in two 
aspects: the comparison of road area extraction result and the comparison of the road centerline extraction 
result. Completeness (COM), Correctness (COR) and Quality (Q) [47] are employed as evaluation metrics 
to measure the both experiments as follows. 

COM = 


TP + FN 


( 21 ) 
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TABLE I 

Quantitative comparisons with state-of-the-art methods in two datasets, where the red values marked in bold 

ARE THE BEST AND THE BOLD BLUE VALUES ARE THE SECOND BEST. (BEST VIEWED IN COLOR) 


Methods 

Imagel 

Image2 

COM 

COR 

Q 

COM 

COR 

Q 

Huang 

0.9479 

0.7954 

0.7621 

0.9339 

0.6142 

0.5886 

Miao 

0.7690 

0.9875 

0.7203 

0.2612 

0.4956 

0.2100 

Shi« 

0.8722 

0.9850 

0.8656 

0.8237 

0.9383 

0.7814 

ShP 

0.8760 

0.9846 

0.8642 

0.8596 

0.9428 

0.8171 

Proposed^ 

0.9264 

0.9997 

0.9261 

0.9223 

0.9172 

0.8514 


TABLE II 

Quantitative comparisons with state-of-the-art methods in Data #1, where the red values marked in bold are the 

BEST AND THE BOLD BLUE VALUES ARE THE SECOND BEST. (BEST VIEWED IN COLOR) 


Methods 

Images 

Image4 

Images 

Image6 

COM 

COR 

Q 

COM 

COR 

Q 

COM 

COR 

Q 

COM 

COR 

Q 

Huang 

0.976 

0.636 

0.627 

0.943 

0.802 

0.765 

0.906 

0.783 

0.724 

0.973 

0.571 

0.562 

Miao 

0.906 

0.780 

0.722 

0.965 

0.853 

0.828 

0.939 

0.913 

0.862 

0.753 

0.601 

0.502 

Shi^" 

0.960 

0.873 

0.843 

0.948 

0.631 

0.610 

0.985 

0.820 

0.810 

0.983 

0.791 

0.780 

Shi^ 

0.920 

0.979 

0.903 

0.871 

0.967 

0.846 

0.871 

0.971 

0.849 

0.891 

0.968 

0.865 

Proposed^ 

0.969 

0.921 

0.894 

0.979 

0.913 

0.895 

0.942 

0.950 

0.897 

0.976 

0.937 

0.917 

Proposed^ 

0.965 

0.937 

0.906 

0.974 

0.917 

0.896 

0.957 

0.950 

0.910 

0.989 

0.944 

0.934 


COR = 


Q = 


TP 


TP + FP’ 
TP 


TP + FN + FP 


( 22 ) 

(23) 


where TP, FP and FN are the true positive, false positive and false negative, respectively. 

For the evaluation of road area extraction result, we compare the extracted road map with reference 
map in the corresponding locations. 

Due to the deviation between the manually labeled centerline and the real centerline, the comparison 
of centerline result is carried out by matching the extracted centerline map to the reference map using 
the so-called “buffer method”, in which every proportion of the network within a given buffer width p 
from the other is considered as matched [48]. That is a predicted centerline point is considered to be a 
true positive if it is at most p distant from a reference centerline point. 


D. Parameter setting 

There are mainly two parameters in the proposed method, i.e. the balancing factor a in GC and the 
scale facor a in TV. In the experiments, we set other parameters as follows: The three different number 
of superpixels are set as 8000, 10000 and 12000 in Proposed^, while in Proposed^ we choose the number 
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(a) Original Image (b) Huang (c) Miao (d) Shi“ ( 6 ) 8111 * (f) Proposed' (g) Proposed* (h) Reference map 



(a) Original Image (b) Huang (e) Miao (d) Shi“ (e) Shi* (f) Proposed' (g) Proposed* (h) Reference map 




Fig. 8. Visual comparisons of centerline extraction results on Data #1. There are four rows and eight columns of subfigures. Each row 
shows the results of one image. From left to right: (a) original image, (b) result of Huang [12], (c) result of Miao [16], (d) result of Shi“ 
[14], (e) result of Shi** [15], (f) result of Proposed^, (g) result of Proposed®, (b) the reference map. 


of superpixels as 10000 in the case of no other special instructions. In MJCR, we randomly choose 50 
positive samples and negative samples from the ground truth. To obtain a sparse coefficients a*, we set 
the weighting factor A = 10. We set the window size of neighborhood as 2cr in NMS. In the fitting based 
connection method, the width of searching window R is set to be a. For Data #1, the road width is 
about 12-15 pixels, thus we choose 3 as the buffer width. While for Data ^2 we set the buffer width p as 
2 pixels. Actually, in the experiments, we find that satisfactory results can be obtained for most images 
when setting a = 0.8 in GC. Thus we keep this parameter fixed in the following experiments. For the 
scale factor a in TV, we get the optimal value for each image via cross validation, which will be discussed 
later. 
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(a) Original Image (b) Huang 


(c) Miao 


(d) Shi" 


(e) Shi' 


(i) Proposed^ (g) Proposed^ (h) Reference map 


(a) Original Image (b) Huang (c) Miao (d) Shi" (e) Shi' (i) Proposed^ (g) Proposed^ (h) Reference map 




Fig. 9. Visual comparisons of centerline extraction results on Data #2. There are eight rows and four columns of subfigures. Each row 
shows the results of one image. From left to right: (a) original image, (b) result of Huang [12], (c) result of Miao [16], (d) result of Shi“ 
[14], (e) result of Shi^ [15], (f) result of Proposed^, (g) result of Proposed^, (b) the reference map. 


E. Comparison of road area extraction results 

In order to illustrate an intuitive comparison of different methods on road area extraction, we display 
the results by comparing the extracted map with the reference map in Fig. 7. Shi“ and Shi^ achieve 
better performance than other two comparing methods, while both of them are inefficient to extract the 
road under the occlusions of cars and trees. Huang’s method can extract the complete road network well. 
However, it brings in more false positives (the red areas in Fig. 7(b)). For Image2, Miao’s method get 
unsatisfactory result. This is because the homogeneous regions are considered to be road after the edge¬ 
filtering method [49], meanwhile in Image2 some forest areas are also homogeneous regions too. Thus 
it is hard to select an appropriate threshold to distinguish them from the road areas. From the comparing 
results, it shows that the proposed method extracted more satisfactory and coherence road map, especially 
for Image2. Besides, our approach is more robust against the occlusions of trees and cars. 

The quantitative performances of road area extraction results are summarized in Tables I. As we know, 
COM and COR can be one-sided metrics. We can enlarge one at the cost of the other one. The Q term 
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combines both of them, thus is a overall measurement. As can be seen from the table, our method achieve 
the highest value in Q. Shi“ and Shi^ obtain higher values than Huang’s method and Miao’s method by 
a large margin. In summary, the Q value of our method is 4% higher than second best method, which 
demonstrates the validity of our approach on the road area extraction. 

F. Comparison of road centerline result on Data #1 

We carried out experiments on all the 30 images in Data #1. Due to space limit, we select 4 of them to 
illustrate the comparing results among the aforementioned methods in Fig. 8. By comparing the extracted 
road centerline map with its corresponding reference map, we can have following conclusions: (1) Through 
morphological thinning based road centerline method (e.g. Huang’s method and Ship’s method) is fast 
and easy to extract road centerlines, this method produces small spurs, which reduce the smoothness 
and correctness of road network. (2) The regression based road centerline methods (e.g. Miao’s method 
and Ship’s method) can extract relatively smooth centerlines well. However, these methods can’t link the 
centerlines well in the intersection areas (see subfigure (c) and (e) of Imagefi in Fig. 8). (3) There are 
some false positives in Miao’s method, this is because some forest areas and bare soil are also tend to be 
homogenous regions in VHR remote sensing images, thus it is hard to choose an appropriate threshold 
to distinguish them with road areas. (4) Our proposed methods (Proposed^ and Proposed^) achieve more 
smooth and continuous road centerlines than other comparing methods, and both of them work well 
in the junction areas. Besides, the centerline result of Proposed^ is more similar to the reference map 
than Proposed^, which demonstrates that multi-scale information fusion is suitable for the road centerline 
extraction task. 

Quantitative comparisons of road centerline extraction result among different methods are illustrated 
in Table II. As we shall see, in terms of COR, the proposed method (Proposed^) achieves the second 
best results after Ship’s method, while in terms of COM, the proposed method obtains bigger values than 
Ship’s method by a large margin. Thus our proposed method achieves relatively higher Q values, which 
is an overall evaluation index. The thinning based methods (i.e. Huang’s method and Ship’s method) and 
Miao’s method give relatively low values of COR and Q due to the appearance of false positives and 
small spurs. Proposed^ obtains the second best performance among all the methods, while it is inferior 
to the Proposed^ almost in each comparing items. 

G. Comparison of road centerline result on Data #2 

Visual comparisons of centerline extraction results on Data #2 among different methods are displayed 
in Fig. 9. we find that there are large areas of dense forests and grasslands in those images, which is hard to 
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Fig. 10. Quantitative comparisons of centerline extraction results with different scale parameters. COM, COR and Q are calculated when 
the buffer with p = 2. (Best viewed in color) 


extract correct road network from them for Miao’s method. Thus the centerline result of Miao’s method 
is incomplete (see Fig. 9(c)). Huang’s method and Ship’s method can extract relatively complete road 
network, while their methods produce small spurs due to the algorithm of centerline extraction. Ship’s 
method and Proposed^ achieve better performance than other methods (e.g. Huang’s method, Miao’s 
method and Ship’s method), while both of them are inferior to the Proposed^ in term of completeness 
and smoothness. 

Table III presents the quantitative comparisons on Data ^2 for different methods. As we shall see, 
Miao’s method achieves low values both in terms of COM and Q, because the road network is hard 
to extract under the conditions of large areas of forests and grasslands. Huang’s method and Ship’s 
method gain relatively large values in COM and small values in COR, which leads to small Q values. 
Ship’s method and Proposed^ achieve similar performance in average among all the four images, while 
Proposed^ achieves higher values than both of them in terms of all the three metrics (e.g. COM, COR 
and Q). Specifically, the average Q value among all the images of the Proposed^ is 3% higher than the 
second best method. 

H. Impact of parameter a 

In the proposed method, tensor voting is used to extract road centerlines. There is only one main 
parameter cr to be tuned in TV. In addition, the window size in NMS and the width of searching window 
in fitting based connection algorithm are also determined by cr. Thus it is important to choose an appropriate 
cr for each image. To give an instructions for users to choose a suitable cr, the sensitivity of o is tested 
in this section. 
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TABLE III 

Quantitative comparisons with state-of-the-art methods in Data #2, where the red values marked in bold are the 

BEST AND THE BOLD BLUE VALUES ARE THE SECOND BEST. (BEST VIEWED IN COLOR) 


Methods 

ImageT 

Images 

Image9 

Image 10 

COM 

COR 

Q 

COM 

COR 

Q 

COM 

COR 

Q 

COM 

COR 

Q 

Huang 

0.804 

0.773 

0.651 

0.872 

0.785 

0.704 

0.893 

0.820 

0.747 

0.941 

0.746 

0.712 

Miao 

0.681 

0.857 

0.612 

0.681 

0.902 

0.634 

0.633 

0.937 

0.607 

0.626 

0.876 

0.575 

Shi“ 

0.889 

0.830 

0.752 

0.911 

0.802 

0.744 

0.864 

0.804 

0.714 

0.939 

0.774 

0.737 

Shi^ 

0.922 

0.977 

0.902 

0.857 

0.972 

0.836 

0.839 

0.944 

0.800 

0.844 

0.940 

0.801 

Proposed^ 

0.927 

0.930 

0.866 

0.950 

0.941 

0.897 

0.905 

0.915 

0.835 

0.949 

0.820 

0.790 

Proposed^ 

0.964 

0.955 

0.922 

0.980 

0.950 

0.924 

0.909 

0.963 

0.878 

0.924 

0.963 

0.892 



Fig. 11. Visual comparisons with different scale parameters for the proposed method. There are two rows and six columns of subfigures. 
From left to right and up to down: (a) original image, (b) classification result of the Proposed^, (c) centerline result with a = 6, (d) centerline 
result with a = 9, (e) centerline result with a = 12, (f) centerline result with a — 15, (g) centerline result with a = 18, (h) centerline result 
with cr = 21, (i) centerline result with a = 24, (j) centerline result with a = 27, (k) centerline result with a = 30, (1) centerline result with 
cr = 33. 

Visual comparisons with different a values are illustrated in Fig. 11. As is seen from the figure, there 
are some diseontinuous regions in the classifieation result (see Fig. 11(b)) due to the occlusions of trees 
(see Fig. 11(a)). One reason to choose TV as centerline extraction algorithm is that TV can eliminate the 
road discontinuity with a suitable scale parameter. When a is small, TV is not only ineffective to link the 
diseontinuous areas, but also ean’t extract centerlines well for those eontinuous regions (see Figs. 11(c), 
11(d) and 11(e)). With the increase of scale parameter a within a certain range, TV can extract centerlines 
for the connected areas and link the centerlines for the disconnected areas (see Figs. 11(f), 11(g), 11(h), 
ll(i) and ll(j)). When a is greater than a eertain value, although TV ean eonneet the diseontinuous 
centerlines well, it also brings in some false positives (see Figs. ll(k) and 11(1)). Those false positives 
lead to wrong navigation information in practical application. Therefore, to get a smooth and complete 
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road network, a suitable value for scale parameter should be carefully selected. 

To give an intuitive comparison, the quantitative performances with different scale parameters are 
summarized in Fig. 10. In this figure, COM, COR and Q are calculated with the buffer width p = 2. 
When a is less than 12 pixels, COM, COR and Q are all increase along with the growth of a. Then with the 
increase of a, COR declines gradually, the overall metric (Q) reduces even though COM remains stable. 
We find that satisfactory results can be obtained with a in a large range (e.g. 15-24 pixels). Specifically, 
the highest Q values are obtained when a is set to be 15 or 18 pixels. It should be noted that the road 
width w of the test image is about 12 pixels, thus it is suitable to set the scale parameter as a = l.bw. 

V. Conclusions 

In this paper, we propose an accurate road centerline extraction method for VHR remote sensing images. 
The proposed method contains three stages: homogeneous road area segmentation, smooth and correct 
road centerline extraction, and centerline connection around the road intersections. Specifically, in the 1®‘ 
stage, to obtain an accurate road area extraction result, multiscale joint collaborative representation and 
graph cuts are introduced to incorporating the mutiscale features and spatial information. In the 2"‘* stage, 
to overcome the shortcomings of morphological thinning algorithm, tensor voting and non-maximum 
suppression algorithm are utilized to extract smooth and correct road centerlines. Finally, in the stage, 
to tackle the ineffectiveness of the existing methods around the road intersections, a fitting based centerline 
connection algorithm is proposed to complete the road network. In terms of both quantitative and visual 
performances, the proposed method achieves better results than all the other comparing methods. Moreover, 
the proposed method are not sensitive to the parameters, which are tuned easily. As another contribution, 
a new and challenging road centerline extraction dataset for VHR remote sensing images will be publicly 
available for further studies. 
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